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Sensor-based degradation signals measure the accumulation of 
damage of an engineering system using sensor technology. Degrada- 
tion signals can be used to estimate, for example, the distribution of 
00 ' the remaining life of partially degraded systems and/or their compo- 

s ! ' nents. In this paper we present a nonparametric degradation model- 

ing framework for making inference on the evolution of degradation 
signals that are observed sparsely or over short intervals of times. Fur- 
thermore, an empirical Bayes approach is used to update the stochas- 
^s . tic parameters of the degradation model in real-time using training 

degradation signals for online monitoring of components operating 
in the field. The primary application of this Bayesian framework is 
CZ2 . updating the residual lifetime up to a degradation threshold of par- 

tially degraded components. We validate our degradation modeling 
approach using a real-world crack growth data set as well as a case 
^. ■ study of simulated degradation signals. 

(N 

r | 1. Introduction. Most failures of engineering systems result from a grad- 

ly-) ■ ual and irreversible accumulation of damage that occurs during a system's 

-^ \ life cycle. This process is known as degradation [Bogdanoff and Kozin (1985)]. 

^^ ■ In many applications, it can be very difficult to assess and observe physical 

degradation, especially when real-time observations are required. However, 
degradation processes are almost always associated with some manifesta- 
tions that are much easier to observe and monitor overtime. Generally, the 
evolution of these manifestations can be monitored using sensor technol- 
ogy through a process known as Condition Monitoring (CM). The observed 
condition-based signals are known as degradation signals [Nelson (1990)] 
and are usually correlated with the underlying physical degradation process. 
Some examples of degradation signals include vibration signals for monitor- 
ing excessive wear in rotating machinery, acoustic emissions for monitoring 
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crack propagation, temperature changes and oil debris for monitoring engine 
lubrication and many others. 

Degradation modeling attempts to characterize the evolution of degra- 
dation signals. There is a significant number of research works that have 
focused on degradation models; these include models presented in Lu and 
Meeker (1993), Padgett and Tomlinson (2004), Gebraeel et al. (2005), Miiller 
and Zhang (2005), Gebraeel (2006) and Park and Padgett (2006). Many of 
these models rely on a representative sample of complete degradation sig- 
nals. A complete degradation signal is a continuously observed signal that 
captures the degradation process of a component from a brand "New State" 
to a completely "Failed State." 

Unfortunately, building a database of complete degradation signals can 
be very expensive and time consuming in applications, such as monitor- 
ing of jet engines, turbines, power generating units, structures and bridges 
and many others. For example, in applications consisting of relatively static 
structures such as bridges, degradation usually takes place very slowly (sev- 
eral tens of years). Since the system is relatively static, it suffices to observe 
the degradation process at intermittent discrete time points. The result is 
a sparsely observed degradation signal such as the signals depicted in Fig- 
ure 1(a). In contrast, in applications consisting of dynamic systems, such 
as turbines, generators and machines, degradation cannot be reasonably as- 
sessed by sparse measurements. At the same time, continuous observations 
of the complete degradation process of such systems are economically unjus- 
tifiable. Usually, the only way to gain a relatively accurate understanding of 
the health/performance of a dynamic system is to monitor its performance 
over a time interval. In naval maritime applications, power generating units 
of an aircraft are removed, tested for a short period of time (during which 
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degradation data can be acquired) and put back into operation. The result 
is a collection of fragmented degradation signals as depicted in Figure 1 (b) . 

In this paper we develop a degradation model that applies to incomplete 
degradation signals as well as complete degradation signals. An incomplete 
degradation signal is defined as a signal that consists of sparse observations 
of the degradation process or continuous observations made over short time 
intervals (fragments). One challenge in such applications is that the evolu- 
tion of the degradation signals cannot be readily assessed to determine the 
parametric form of the underlying degradation model. This is because one 
cannot clearly trace how a degradation signal progresses over time from in- 
complete observations. For example, is there a well defined parametric model 
that describes the signals in Figure 1? 

To overcome this challenge, the underlying degradation model in this pa- 
per is assumed nonparametric. Most degradation models used to character- 
ize the evolution of sensor-based degradation signals are parametric models. 
A common approach is to model the degradation signals using a parametric 
(linear) model with random coefficients [Lu and Meeker (1993); Gebraeel 
et al. (2005); Gebraeel (2006)]. Other modeling approaches assume that the 
degradation signal follows a Brownian motion process [Doksum and Hoyland 
(1992); Pettit and Young (1999)] or a Gaussian process with known covari- 
ance structure [Padgett and Tomlinson (2004); Park and Padgett (2006)]. 
In contrast, we assume that the mean and covariance functions of the degra- 
dation process are unknown and they are estimated based on an assembly 
of training incomplete degradation signals. The mean function is estimated 
using standard nonparametric regression methods such as local smoothing 
[Fan and Yao (2003)]. The covariance function is decomposed using the 
Karhunen-Loeve decomposition [Karhunen (1947); Loeve (1945)] and esti- 
mated using the Functional Principal Component Analysis (FPCA) method 
introduced by Yao, Miiller and Wang (2005). 

Under the nonparametric modeling framework, one condition for accurate 
estimation of the mean and covariance functions is that the degradation pro- 
cess is densely observed throughout its support. However, in applications 
where the degradation signals are incompletely sampled, not all degrada- 
tion signals are observed up to the point of failure; in addition, only a few 
components will survive up to the maximum time point of the degradation 
process support. Consequently, the degradation process is commonly under- 
sampled close to the upper bound of its support. To overcome this difficulty, 
we introduce a nonuniform sampling procedure for collecting incomplete 
degradation signals, which ensures relatively dense coverage throughout the 
sampling time domain. 

One important application of degradation modeling is predicting the life- 
time of components operating in the field. For real-time monitoring, an em- 
pirical Bayes approach is introduced to update the stochastic parameters of 
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the degradation model. In this paper we focus on estimation of the distribu- 
tion of the residual life up to a degradation threshold for partially degraded 
components using training degradation signals which are sparsely or com- 
pletely observed. Other applications of the degradation modeling and the 
Bayesian updating procedure are estimation of the lifetime at a specified 
degradation level and estimation of the degradation level at a specified life- 
time. 

We evaluate the performance of our methodology using both a crack 
growth data set and simulated degradation signals. In these empirical studies 
we compare parametric to nonparametric degradation modeling, assess the 
estimation accuracy of the remaining lifetime for complete and incomplete 
signals, and contrast uniform vs. nonuniform sampling procedures for acquir- 
ing ensembles of incomplete degradation signals. In both studies there is not 
a significant decrease in the accuracy of the residual life estimation when 
using ensembles of incomplete instead of complete signals. We also highlight 
the robustness of our approach by comparing it with misspecified paramet- 
ric models, which are common when the underlying degradation process is 
complicated and sparsely observed. Last, we show in the simulation study 
that using a nonuniform sampling procedure that ensures dense observation 
of the sampling time domain reduces the estimation error. Based on these 
empirical studies, we conclude that the nonparametric approach introduced 
in this paper is efficient in characterizing the underlying degradation process 
and it is more robust to model misspecification than parametric approaches, 
which is particularly important when the training signals are incompletely 
observed (sparse or fragmented). 

The remainder of the paper is organized as follows. Section 2 discusses the 
development of our degradation modeling framework. The empirical Bayes 
approach for updating the degradation distribution of a partially degraded 
component is introduced in Section 3. The derivation of the remaining life- 
time distribution under the empirical Bayes approach is presented in Sec- 
tion 4. In Section 5 we introduce an experimental design for sampling incom- 
plete degradation signals. We discuss performance results of our methodol- 
ogy using real-world and simulated degradation signals in Section 6 and 7, 
respectively. 

2. Sensor-based degradation modeling. We denote the observed degra- 
dation signals Sj(ijj), for j = l,...,mj (m; is the number of observation 
time points for signal i) and i = 1, . . . , n (n is the number of signals) where 
{tij}j=i,...,mi are the observation time points in a bounded time domain 
[0,M] for signal i. Note that M will always be finite since any industrial 
application has a finite time of failure. We model the distribution of the 
signals Si(t) by borrowing information across multiple degradation signals. 
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We decompose the degradation signal as 

(2.1) S i (t)=fi(t) + X i {t) + ae l {t), 

where fj,(t) is the underlying trend of the degradation process and is as- 
sumed to be fixed but unknown, and Xi(t) represents the random deviation 
from the underlying degradation trend. We also assume Xi{t) and £«(£) are 
independent. 

The model in (2.1) is a general decomposition for functional data with 
various modeling alternatives and assumptions for the model components: 
/i(t), Xi{t), and £i(t). In this paper we discuss one such modeling alterna- 
tive which applies to sparse and fragmented signals as well as to complete 
signals and it applies under the assumption that the observation time points 
{tij}j=i,...,mi are fixed but not necessarily equally spaced and the assump- 
tion that the error terms £i{t) are independent and identically distributed. 
Deviations from these assumptions may require some modifications to the 
modeling approach discussed in this paper. 

In our modeling approach, the degradation signal Si(t) follows a stochas- 
tic process with mean /j,(t) and stochastic deviations Xi{t) with mean zero 
and covariance cov(i,t'). The mean function /j,(t) and the covariance surface 
cov(t,t') are both assumed to be nonparametric, that is, no prespecified 
assumption on their shape. This generalized assumption encompasses the 
particular cases developed earlier by Gebraeel et al. (2005) and Gebraeel 
(2006), which assume a linear trend, fi(t) = a + j3t where a ~ iV(0, 8 a ) and 
/3 ~ iV(0, (5^), and parametric covariance structure cov(i, t') = 6 a + Sptt' '. 

The following steps discuss how we estimate the mean function and the 
covariance surface of our degradation model. 

Step 1: We use nonparametric methods to estimate the mean /u(i). In 
this paper we use local quadratic smoothing [Fan and Yao (2003)] to allow 
estimation of the mean function under general settings including complete 
and incomplete (sparse and fragmented) signals. The bandwidth parameter, 
which controls the smoothing level, is selected using the leave-one-curve-out 
cross-validation method [Rice and Silverman (1991)]. Alternative estima- 
tion methods include decomposition of the mean function using a basis of 
functions (e.g., splines, Fourier, wavelets) and estimate the coefficients us- 
ing parametric methods. These alternative methods will apply under various 
signal behaviors (e.g., smooth vs. with sharp changes, uniformly vs. nonuni- 
formly sampled). 

Step 2: The covariance surface is estimated using the demeaned data, 
Si(t) — fi(t), where fi(t) is the local quadratic smoothing estimate of /u(t). 
Using the Karhunen-Loeve decomposition [Karhunen (1947); Loeve (1945)], 
the covariance, cov(t,t') = Cov(Sj(i) , Si(t')) , can be expressed as follows: 

oo 

(2.2) cov(t,t') = J2*kfa(t)Mt'), t,t'e[0,M], 

fe=i 



(2.4) XiiUj) =J2&kMti 
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where (f>k(t) for k = 1,2, ... are the associated eigenfunctions with support 
[0, M] and Ai > A2 > • ■ • are the ordered eigenvalues. Based on this decom- 
position, the deviations from the underlying degradation trend Xi(t) are 
decomposed using the following expression: 

00 
(2.3) Xi(tij) =y^^ik(f>k(tij), 

fc=l 

where ^ called scores are uncorrelated random effects with mean zero and 
variance E(£? fc ) = A&. The decomposition in equation (2.3) is an infinite sum. 
Generally, only a small number of eigenvalues are commonly significantly 
nonzero. For the eigenvalues which are approximately zero, the correspond- 
ing scores will also be approximately zero. Consequently, we use a truncated 
version of this decomposition. Therefore, expression (2.3) can be approxi- 
mated as follows: 

K 

h i] ) 1 
fe=l 

where K is the number of significantly nonzero eigenvalues. We select K to 
minimize the modified Akaike criterion defined by Yao, Miiller and Wang 
(2005). 

In the statistical literature this method has been coined Functional Prin- 
cipal Component Analysis (FPCA). The key reference for FPCA is Ram- 
say and Silverman [(1997), Chapter 8]. Another important reference is Yao, 
Miiller and Wang (2005), in which the authors derived theoretical results for 
model parameter consistency and asymptotic (n large) distribution results 
under the assumption that the scores follow a normal distribution. 

An alternative method for estimating the covariance function of the pro- 
cess Xi(t) is decomposing the covariance function as in equation (2.2) where 
the basis of functions {(pk,k = 1,2,...} is fixed [James, Hastie and Sugar 
(2000)]. However, this approach doesn't allow dimensionality reduction in 
the same way FPCA does and it is not theoretically founded. 

3. Degradation model updating. Next, we consider a component oper- 
ating in the field called fielded component. Assume that we have observed its 
degradation signal at a vector of time t = (£1, . . . ,t m *); therefore, S(t) de- 
notes the observed signal of the testing component, m* represents the num- 
ber of observations and t* = t m * denotes the latest observation time. In this 
section we introduce an Empirical Bayes approach which allows real-time up- 
dating of the distribution of the degradation process for partially degraded 
components given the observed signal S'(t) and the prior distribution of the 

scores £jj. for k = 1,2, The prior distribution of the scores is estimated 

empirically from a set of historical degradation signals. 
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Proposition 1 illustrates the updating procedure assuming that the prior 
distribution of the scores is normal and assuming that the mean function /i(t) 
and the basis of functions (f>k(t), k = 1, . . . ,K, are fixed. The proof of this 
proposition follows from the theory of Bayesian linear models. 

Proposition 1. Assume that S(t) follows 

K 



S(t) = fi(t) +J2^k(-t) + <t), 



k=l 

where the prior distribution of £& is N(0,Xk) with £i,...,£k uncorrelated; 
e(t) are independent of^k for k = 1, . . . , K; the distribution of s(t) is N(0, a 2 ) 
with a 2 fixed. It follows that the posterior distribution of the scores is 

(&...,&)' ~N(Cd,C), 

where 

C=(-^P(t)'P{t) + A- 1 } and rf=i?(t)'(S(t)-^(t)) 

with 

S(t) = (S(h), ..., S(t m *)Y, M (t) = (m(*i), • • • ,M*m.))', 
(3-1) 

/ Mh) ■■■ <t>K(h) 
A = diag(Ai,...,A^), P(t)= ... 

In Proposition 1 the prior distribution of the scores is specified by the 
variance parameters A&, k = 1,. . . ,K, which are estimated using the degra- 
dation model in Section 2 and based on a set of incomplete or complete 
training degradation signals. Specifically, we first apply Functional Princi- 
pal Component Analysis on the historical degradation signals which will 
further provide estimates for the variance parameters A^, k = l, . . . ,K, and 
the eigenfunctions (ftk, k = 1, . . . , K. Based on these estimates, we obtain the 
posterior distributions of the updated scores ££,. . . ,£,* K since the matrix C 
and the vector d are fully determined by the eigenvalues Afc, k = 1, . . . ,K, 
and the eigenfunctions (pk, k = 1, . . . ,K. The expectation of the posterior 
scores is nonzero and, therefore, we denote the posterior mean function 

Following Proposition 1, the expectation of the posterior distribution fol- 
lows the same formula as the conditional expectation estimator in Yao, 
Miiller and Wang (2005), equation (4). Generally, this similarity applies 
under the empirical Bayesian prior derived from FPCA. On the other hand, 
the sampling distribution of the conditional expectation estimator in Yao, 
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Miiller and Wang (2005) is different from the posterior distribution of ££, k = 
1,...,K, since their variances are not equal. Moreover, the conditional ex- 
pectation estimator and its mean estimation error in Yao, Miiller and Wang 
(2005) are conditional on the training observations, whereas the posterior 
distribution in Proposition 1 is conditional on the observations of a new 
component. 

The advantage of this Bayesian framework is that it unifies the condi- 
tional expectation estimation and prediction into a procedure which allows 
updating the distribution of the degradation process for a new component. 
We can therefore use the posterior distribution of the scores for a partially 
degraded signal to estimate the distribution of various statistical summaries, 
including the lifetime at a specified degradation level and estimation of the 
degradation level at a specified time. In the next section we discuss one 
specific application to this updating framework: residual life estimation. 

4. Remaining life distribution. In this paper we focus our attention on 
engineering applications where a soft-failure of a system occurs once its 
underlying degradation process reaches a predetermined critical threshold. 
This critical threshold is commonly used to initiate maintenance activities 
such as repair and/or component replacement well in advance of catastrophic 
failure. Consequently, degradation data can still be observed beyond the 
critical threshold. In this section we describe how our degradation modeling 
framework is applied to estimate the distribution of remaining life up to 
a degradation threshold of partially degraded systems. 

In the remainder of this section, S*(-) will denote the underlying degra- 
dation process of a partially degraded system. Based on the degradation 
process S*(-), the failure time of a system is defined as 

(4.1) T= inf ■ {S*(t)>D}. 

te{o,M] 

One has to bear in mind that T may not exist if the threshold D is set too 
high, that is, the component may fail before its degradation signal reaches 
the threshold. The selection of the failure threshold D is an important prob- 
lem, however, this aspect is beyond the scope of this paper. In this work, we 
assume that T exists, and the threshold D is known a priori. This is a rea- 
sonable assumption because in many industrial applications failure/alarm 
thresholds are usually based on subjective engineering judgement or well- 
accepted standards, such as the International Standards Organization (ISO) 
(e.g., the ISO 2372 is used for defining acceptable vibration threshold levels 
for different machine classifications) . A second assumption is that the failure 
time T is smaller than a maximum failure time M. This assumption is also 
reasonable, as in practice a component may be replaced after a given period 
of time even if it did not fail. 
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The distribution of the residual life (RLD) of a partially degraded compo- 
nent at a fixed time t* S [0, M] is estimated assuming that the degradation 
process S*(-) of the component follows a posterior distribution based on 
Proposition 1. We estimate the distribution of the residual life (RLD) using 

R(y\t*) = P(T - t* < y\S*(t) ~ Gaussian^* (£), Cov*(i, t')), t* <T<M), 

where //*(£) and Cov* (£,£') are the posterior mean and covariance functions 
of the degradation process £*(•). The derivations of /x*(t) and Cov* (£,£') are 
based on the results of Proposition 1. We note here that the distribution 
of the RLD above is not conditional on the observed signal of the partially 
degraded component but on the posterior distribution of its degradation 
process; since the degradation is only partially observed and most often 
sparsely sampled, conditioning on the posterior distribution will generally 
provide a more accurate RLD estimator since we incorporate the additional 
information in the training degradation signals. 

Furthermore, we estimate RLD under two assumptions: 

A.l The new component has not failed up to the last observation time point 
£*, that is, the failure time becomes 

T= inf {S*(t)>D}= inf \S*(t) > D\ := T*. 
te[o,M] te[t*,M] 

A.2 We assume the probability that the degradation process S*(t) crosses 
back the threshold D after the failure time T* is negligible, that is, 
P(S*(T* + y) < D) ~ for all y > 0. This implies, if we condition on 
y > T* — t* > 0, which is the same as conditioning on T* < t* + y, 
P(S*{t* +y)< D\T* < t* + y) w 0. This further implies 

P(S*{t* + y)>D) = P{S*{t* + y)> D\T* < t* + y)P{T* <t* + y) 

&P(T*<t* + y). 

Under these two assumptions, the RLD becomes 

R(y\t*) = (by A.l) P(T*-t* < y\S*(t)) « (by A.2) P(S*{t* +y)> D\S*(t)). 

The approximation in assumption A.2 is similar to the approximation in 
the paper by Lu and Meeker (1993) which assumes that the probability of 
a negative random slope in the linear model is negligible. One particular 
case for the assumption A.2 to hold is that the signal is monotone. However, 
monotonicity is not a necessary condition. Assumption A.2 also holds for 
nonmonotone signals — an example of such a signal is in Figure 2. 

Proposition 2 below describes the updating procedure for RLD of a new 
component given the posterior distribution of its degradation process S*(-) 
updated up to time t*. The proof follows directly as a consequence of Propo- 
sition 1. 
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Observed 




t * T * t*+y 

Failure time 

Fig. 2. Example of a signal for which assumption A. 2 holds. 



Proposition 2. For a new partially degraded component with its degra- 
dation process S*(-) updated up to time t* , the residual life distribution is 
given as follows: 



(4.2) P(T-t* <y\S*(-),T>t* 



$z(g*(y\t*))-$z(g*(o\t*)) 



l-*zO/*(0|t*)) 

where &z represents the standard normal cumulative distribution function 
and g*(y\t*) = »'f + y^ D with 

fi*(t* +y) = n(t* + y) + (Cd)'p{t* + y), 

K K 

V*(t* + y)=^2 5D[C*i,k2fox(<* +V)<t>k 2 (t* +y)]. 
fe 1= ife 2 =i 

In the above equations, p(t* + y) = {(/>i(t* + y), . . . , <f>K(t* + y))', and Cki t k2 
refers to the (^1,^2) element of the matrix C. 

One advantage of obtaining the distribution rather than simply a point 
estimate is that we can also derive a confidence interval for the remain- 
ing lifetime up to a degradation threshold D. Following the derivation in 
Proposition 2, a 1 — a confidence interval for RLD is [L, U] such that 

P(L <T-t* < U\S*(-),T >t*) = l-a. 
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Since we have one equation with two unknowns, the lower — L and the 
upper — U tails are commonly equally weighted, and, therefore, 

<S> z (g*(U\t*))-<S> z (g*(0\t*)) =l a 
l-$z(s*(0|t*)) " 2 

and 

$z(g*(L\t*))-$ z {9*{0\t*)) = a 
l-$z(<7*(0|<*)) 2' 

However, we cannot obtain exact solutions for L and U because we do not 
have a closed-form expression for the inverse of the cumulative density func- 
tion of T — t* . For example, the first relationship is equivalent to finding U 
from g*(U\t*) = z ai where z ai is the 1 — a\ quantile of the normal distribu- 
tion. Using this equation, we would like to obtain U such that 

fi*(t* + U)-D _ 

^/v^WtuT ~ Zai ' 

which is a nonlinear function of U and its solution does not have a close form 
expression. We therefore resort to parametric bootstrap [Efron and Tibshi- 
rani (1993); Davison and Hinkley (1997)] to sample from the distribution 
of T — t* which will give us a set of realizations from this distribution — 
Ti,T2, . . . ,Tb- Using these realizations from the distribution of T — t*, we 
estimate a quantile bootstrap confidence interval. 

The confidence interval estimation procedure is as follows. For b = 1, . . . , B: 

1. Sample £ = (£^, . . . , ^ K ) from the multivariate normal distribution of the 
posterior scores provided in Proposition 1. 

2. Obtain a simulated signal 

K 



s b (t)=v(t)+j2&<t>k(t), 



fc=i 

where £|, k = 1,. . . ,K, are the scores sampled at Step 1. 
3. Take T b = mf te[0M] {S b (t) > D}. 

Using the sampled values T\ , T2, . . . , Tb , we compute the empirical a/2 and 
(1 — a/2) quantiles, T a j 2 and T 7 1 _ a / 2 , respectively. We estimate the upper 
and lower bound of the confidence interval by L = T a / 2 and U = T 1 _ a / 2 - It 
follows that [L, U] is an approximate 1 — a quantile bootstrap confidence 
interval for the residual life time of the fielded component. 

An additional approach to the (parametric) bootstrap method described 
above is to (re)sample the signal data resulting in multiple bootstrap sam- 
ples. For each bootstrap sample, estimate the residual lifetime using the 
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approach discussed in this paper; therefore, we obtain a set of realization 
from the distribution of T — t*. In contrast to the bootstrap method de- 
scribed above, this alternative bootstrap approach requires estimating the 
FPCA model for each bootstrap sample which is computationally expensive. 

5. Sampling scheme. The nonparametric degradation modeling frame- 
work introduced in this paper applies to both complete as well as incomplete 
degradation signals. For applications involving incomplete degradation sig- 
nals, it is important to develop a sampling plan that ensures accurate estima- 
tion of the mean function and the covariance surface. Yao, Miiller and Wang 
(2005) provide theoretical results on the estimation of the covariance surface 
using FPCA under large n but small m; for i = 1, . . . , n. In other words, for 
these results to hold, the observation time points {tij}j=i mi i=i n need 
to cover the time domain, [0, M], densely. 

Using the traditional uniform sampling technique, the number of obser- 
vations per time interval decreases as more signals fail, leading to an unbal- 
anced number of observations per time interval — more observations at the 
beginning of the observation time domain but fewer observations at the end 
of the time domain. Further, this unbalanced design will result in decreasing 
estimation accuracy (higher variances) of the mean and covariance estimates 
at later time points. In order to balance the number of observations per time 
interval throughout the time domain [0, M], we propose an experimental de- 
sign using nonuniform sampling. The proposed technique ensures relatively 
dense coverage of the sampling time domain, [0, M] , where M represents the 
last observation time of the longest possible degradation signal for a given 
application. 

We note here that the sampling technique requires input of M at the 
beginning of the experiment although M is unknown. It is often the case 
that, in practice, an experimenter will set a timeline at the beginning of the 
experiment which will specify a limit of how long the experiment will be 
run (e.g., one year vs. one month). This upper limit will specify M. Gen- 
erally, starting with a lower initial value for M will allow the experimenter 
to sample densely enough while having the option to update the sampling 
technique (update M) if not all training signals have reached the failure 
threshold by the initial value for M . 

The following steps outline a sampling procedure for obtaining sparsely 
observed and fragmented degradation signals: 

Step 1: We begin by performing nonuniform sampling of the time domain 
[0, M], thus obtaining a sequence of time points, = t\ < £2 < • ■ • < t m -\ < 
t m = M, for large m. Since only a few components will survive up to the 
maximum time point, M, we increase the sampling frequency at later time 
points in order to cover the sampling time domain at the extreme point, M. 
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Consequently, we sample exponentially, that is, the time interval between 
two consecutive sampling time points decreases exponentially over time (the 
decreasing rate is implicitly determined by the value of M and the number 
of sampling time points). 

Step 2: This step provides a potential sampling timetable (or monitoring/ 
observation schedule) for sparsely observed and fragmented degradation sig- 
nals. We begin by selecting n components. For each component, we select 
its sampling time points from the set t\,...,t m without any prior knowledge 
about their degradation process and lifetime. Next we define two settings: 

Setting 1: This setting is used to obtain sparsely observed degradation 
signals. For component i, we randomly sample rrii time points from the set 
of total time points {t±, . . . , t m }. This results in a set of sparse sampling time 
points, {tn, ■ ■ ■ lUrrn} for this component. 

Setting 2: This setting is used to obtain fragmented degradation signals. 
Recall that fragmented signals are obtained by continuously monitoring 
a component over a short time interval, hence the term "fragment." For 
component i, we select two or more time points B\,B2,... from the set of 
total time points {t\,. . . ,t m }. These points represent the beginning times 
of the signal fragments or sampling intervals. The duration of the sampling 
interval will depend on the type of application, the availability of monitor- 
ing/testing equipment and the associated costs/economics. Consequently, 
the end time points, E\, E2, ■ ■ ■ , will vary from one experiment to another. 
In other words, for component i, we may have two or more time intervals: 
[B^ijE^i], [Bi t 2,Ei t 2\ — 

Step 3: Finally, we observe the degradation signal for the selected com- 
ponents at the selected time points according to the type of incomplete 
signals, sparse or fragmented, and obtain the set of sampled signals Si(tij) 
for i = 1, . . . , n and j = 1, . . . , m*. 

It is important to stress that we select the sampling time points in Step 2 
before observing the degradation signals. Since we do not observe the failure 
time before selecting the time points, we cannot ensure that the degrada- 
tion signal will be observed for all selected sampling time points. This is 
because some components may fail before the latest selected time point. 
Consequently, the observation time points are a subset of the sampling time 
points and will be less densely sampled close to M since the missing observa- 
tions (the difference set between sampling and observation time points) will 
increase in density closer to the upper bound M. In Figure 3 we compare 
the sampling time points selected at Step 2 to the observation time points 
for three components. In this example the sampling time points are nonuni- 
formly selected, whereas the observation time points are approximately uni- 
form since for the first two components, we do not observe at the latest 
times — only the third component fails after its latest sampling time. 
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Two parameters that are used for tuning the sampling plan are as fol- 
lows: the total number of sampling time points m to tai = ]CILi m * an d the 
total number of components n. The more sparsely the signals are observed 
(mtotai is small) , the more signals we need to observe (n needs to be large) . 
Selecting n and mtotai optimally is important to ensure accurate modeling 
of the degradation process at a feasible cost. The larger the number of com- 
ponents n and/or the larger the number of time points mtotai are the higher 
the costs associated with monitoring and testing. Note that selection of n 
and mtotai wm vary according to the type of application. 

6. Case study: Crack growth data. In this section we study crack growth 
data that can be found in various domains of engineering applications, such 
as infrastructure (bridges, steel structures), maritime (hulls of oil tankers), 
aeronautical (aircraft fuselage), energy (vanes of gas turbines), etc. We con- 
sider a situation in which crack growth data can be observed from identical 
units (say, several ship hulls, or turbines) up to a predetermined time pe- 
riod, denoted by M in this paper. A constant threshold, D, is a critical 
crack length representing a soft failure when maintenance and repair should 
be performed. Within this context, we assume that catastrophic failure, that 
is, hard failure, may occur at a relatively larger crack length. 

The data set used in our case study was first published in Virkler, Hill- 
berry and Goel (1979), and has been previously analyzed in other journal 
articles [Kotulski (1998); Cross, Makeev and Armanios (2006) and the ref- 
erences therein]. The specimens in the test were 2.54-mm-thick and 152.4- 
mm-wide center cracked sheets of 2024-T3 aluminum. The crack propagation 
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signals of these specimens were recorded under identical experimental con- 
ditions. In this data set, the crack length was measured in millimeters and 
the observation time was measured by the cumulative load cycles. More de- 
tails about this data set can be found in Virkler, Hillberry and Goel (1979). 
In this study, we set the soft failure threshold to D = 27 mm. We provide 
additional results for another soft threshold in the supplemental material 
[Zhou, Serban and Gebraeel (2010)]. To be consistent with the methodology 
in this study, the observations are censored at common value M = 230,000 
cycles. A representative example of sparsely sampled degradation signals is 
in Figure 4(a). 

6.1. Results and analysis. We report the prediction accuracy of the re- 
maining life for varying time points t* defining the latest observation time 
of a partially degraded component. We consider the following degradation 
percentiles: 10% (the signal has been observed up to time t* , which equals 
to 10% of the lifetime), 20%, . . . , 80% and 90%. For each crack, we pre- 
dict the updated residual lifetime at each of the nine percentiles using the 
degradation signal observed up to that respective percentile. The number 
of signals in this study is 59. We randomly select 50 of the total signals 
as training signals for estimating the model components, and the rest are 
validation signals for evaluating the performance of our model in predicting 
residual life. For each validation signal, we use the following error criteria to 
assess the prediction accuracy: 

I Estimated Life — Actual Lifel 

(6.1) error = — — — — . 

Actual Life 

We replicate the above procedure for 100 times, and report the distribution 
of the errors across the 100 simulations using a set of boxplots, each boxplot 
corresponding to a degradation percentile for the testing components and 
providing the absolute prediction errors for that percentile. 

We first discuss the performance of our nonparametric model for com- 
plete, sparse and fragmented degradation signals. In each complete degra- 
dation signal, we have about 50 observations per signal. To obtain a sparsely 
observed degradation signal, we randomly sample m = 6 observations from 
each complete signal. We use two intervals per signal to obtain fragmented 
degradation signals. The results are illustrated in Figure 4(b)-(d). The re- 
sults indicate that our nonparametric model performs well for complete as 
well as incomplete degradation signals, and the performance is better when 
the incomplete degradation signals are sparse rather than fragmented. Al- 
though we have only approximately 10% observations of complete degrada- 
tion signals under the sparse sampling scenario, the prediction errors do not 
increase significantly. This observation is important in practice; under bud- 
get limitations, one may resort to sparse or fragmented degradation signals 
without significant loss of predictive capability. 
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Fig. 4. The prediction error of residual life prediction for the crack growth data set. 
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We also demonstrate the benefits of our proposed nonparametric degra- 
dation model by comparing it with parametric models as benchmarks. Since 
the degradation signals have a nonlinear trend with a curvature similar to 
the exponential function, we transform the degradation signals using the 
natural logarithm in order to linearize the trend and then apply a linear 
random effects model (henceforth, denoted by "log-linear"). Since under the 
log-transform model, the residual life predictions are inaccurate compared 
to the nonparametric approach, we consider a double logarithm transfor- 
mation of the degradation data (henceforth, denoted by "log-log- linear" ) . 
The results of the sparse scenario using the parametric models "log-linear" 
and "log-log-linear" are reported in Figure 4(e)-(f), respectively. We find 
that both parametric models provide less accurate predictions of the resid- 
ual life than our nonparametric model. This is due to the inaccuracy of the 
parametric models in capturing the crack propagation trend. 

We provide one example in Figure 5 to illustrate the source of the bias 
of the "log-linear" model. In this figure the x-axis represents the degrada- 
tion time and the y-axis represents the crack length, but in the log scale. 
We have one complete, sparse and fragmented degradation signal in Figu- 
re 5(a)-(c), respectively. If the "log-linear" model is the true underlying 
parametric model, we should see a linear trend in all three plots. This seems 
to be true in the sparse or fragmented cases [see Figure 5(b)— (c)]. However, 
for Figure 5(a) showing a complete signal, we note that the degradation trend 
is still nonlinear; the log-transformation does not linearize the signal (the 
same applies for the "log-log" transformation). Therefore, the "log-linear" 
model does not accurately capture the crack propagation trend throughout 
the unit's lifetime. This example shows the potential difficulty of identify- 
ing a reasonable parametric model for sparse and fragmented degradation 
signals and, in turn, demonstrates the robustness of our proposed nonpara- 
metric model to model misspecification. 

7. Simulation study. In this section we simulate nonlinear degradation 
signals from three different models to demonstrate the benefits of using our 
proposed nonparametric degradation modeling approach. We evaluate our 
approach in terms of the prediction accuracy of estimating the residual life 
for complete, sparse and fragmented degradation signals, contrast uniform 
and nonuniform sampling procedures for acquiring the ensembles of incom- 
plete degradation signals, and also investigate the robustness of our model 
to violations of its model assumptions. 

7.1. Simulation models. The degradation signals are simulated from three 
different models, and all of them are special cases of the general model (2.1). 
More specifically: 
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• In Model 1, we choose fx{t) = 30t 2 , Xi(t) = £i<M0, where £i ~ N(0, f ), 
<j )l (t) = \fht 2 , 0<£< 1, and a = l. 

• In Model 2, we choose fi(t) = 30t 2 , Xi(t) = £\<j>\{t) + ^h(t), where £i ~ 
iV(0,3 2 ), <6i(t) = 2t, and £ 2 ~ AT(0, (|) 2 ), <fa(t) = \/80t 2 -|\/80t, < t < 1. 
(The coefficients of the eigenfunctions are chosen so that they form an or- 
thonormal functional basis for < t < 1.) 

• In Model 3, we choose //(*) = 30i 2 - 2sin(47rt), Xi(t) = £i(f>i(t) + 6</>2(»> 
where 6 ~ ^(O^ 2 ), fa(t) = 2t, and £2 ~ 7V(0,(§) 2 ), fo{t) = \/tot 2 - 
|\/80t, < t < 1. 

We simulate from Model 1 because its residual life distribution can be eas- 
ily derived from training signals and updated using validation signals using 
the procedure in Gebraeel et al. (2005). The derived residual life distribu- 
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tion can then be utilized as a benchmark to assess the performance of our 
nonparametric approach. 

Across all the models, the failure threshold is set to D = 10. We gener- 
ate n = 100 "training" signals and n = 100 "validation" signals from each 
model. For a complete signal, we have 51 observations made at an equally 
spaced grid cq, . . . , C50 on [0, 1] with cq = 0, C50 = 1. A sparse or fragmented 
signal is then sampled from a complete signal such that we observe about 6 
observations per signal. The stopping time for each training signal (the last 
point at which a signal is observed) is generated from Uniform distribution 
[Uniform(0.7, 1)] — our simulation results are insensitive to the selection of 
the stopping time distribution. 

We run simulations for 100 times. For each simulation, we compute the 
prediction errors at the following degradation percentiles: 10%, 20%, . . . , 
70%, 80% and 90% of the simulated degradation signals. 

7.1.1. Results and analysis of Model 1. In Figure 6(b)— (d), we present 
the boxplots of the prediction errors when using the nonparametric degra- 
dation model in this paper for complete, fragmented and sparse degradation 
signals. For the sparse scenario, we compare the prediction accuracy of using 
the true parametric model [see Figure 6(e)] and our nonparametric model 
when signals are uniformly sampled [see Figure 6(f)] or nonuniformly sam- 
pled [see Figure 6(d)]. We assess the robustness to model assumptions by 
simulating signals from the model with £1 following a Gamma or Student t 
distribution [see Figure 6(g)-(h)]. We also compute the prediction errors 
under different error distributions [see Figure 6(i)]. 

The first observation is that there is insignificant difference in the predic- 
tion errors between the true parametric model and the nonparametric degra- 
dation model. The differences are larger for high degradation percentiles. 
Since the difference in the prediction errors increases with additional data, 
we observe for a new component, we infer that this small inefficiency arises 
due to a decreased accuracy in the estimation of the empirical prior distri- 
bution at the later time points. 

The second important observation is that the nonuniform sampling tech- 
nique proposed in Section 5 enhances the prediction accuracy of the residual 
life. In Table 1 we list the median prediction errors based on nonuniform 
sampling and uniform sampling techniques. The first row of this table rep- 
resents the time percentile of the degradation signals used for predicting the 
residual life. It is apparent that the nonuniform sampling technique provides 
smaller prediction errors, especially at high time percentiles. This is because 
nonuniform sampling ensures dense coverage of observations over the whole 
time domain, including the region near maximum observation time (M), and 
hence provides more accurate estimate of the mean and covariance functions 
of the model, especially at higher time percentiles. 
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Fig. 6. The prediction error of the residual life estimate for Model 1. 
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Fig. 6. (Continued). 

Last, we assess the robustness to departures from our model assumptions: 
normality of the scores and normality of the errors. In Figure 6(g)— (h), we 
compare the prediction errors when the scores follow Gamma and Student t 
distribution. We also present the results when the errors follow Student t 
distribution in Figure 6(i). The prediction errors for all these different set- 



Table 1 
Prediction errors based on sparse degradation signals that are uniformly or nonuniformly 

sampled 
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Fig. 7. Confidence interval estimation: the coverage rate (a) and mean length (b). In 
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respectively. 



tings are similar. This robustness property of our degradation modeling is 
inherited from the robustness of the FPCA method [Yao, Miiller and Wang 
(2005)]. 

We also evaluate the accuracy of the confidence interval estimates intro- 
duced in Section 4. In Figure 7 we present the coverage rate level and the 
mean of the confidence interval length at the degradation lifetime percentiles 
50%, 60%, 70%, 80% and 90%. The confidence interval level is 1 - a = 0.9. 
The coverage rate is higher for complete signals than for sparse signals 
throughout all percentiles, but the difference is insignificant. The cover- 
age rate for both complete and sparse signals is approximately equal to the 
confidence level 1 — a = 0.9. Moreover, the mean length decreases for higher 
percentiles, implying that the accuracy of the residual life estimate increases 
as the latest observation time point t* is closer to the failure time. 

7.1.2. Results and analysis of Model 2. In the following analysis we still 
use Model 1 as the assumed parametric model and its derived residual life 
distribution as the benchmark. This assumed parametric model correctly 
captures the mean degradation trend of Model 2 but not the underlying 
covariance structure of the degradation process. It is worth mentioning that 
most existing parametric approaches focus on identification of the functional 
form for the underlying degradation trend, ignoring the underlying covari- 
ance structure. 

The results in Figure 8 indicate that our nonparametric model is more 
accurate than the assumed parametric model in predicting the residual life. 
This is because our proposed nonparametric approach, which is FPCA- 
based, cannot only estimate the mean trend accurately but also capture 
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Fig. 8. The prediction error of the residual life estimate for Model 2. 



the dominant modes of the covariance structure correctly. In contrast, para- 
metric models are not flexible enough to accurately capture the underlying 
covariance structure. 

We also compute the prediction error results for cases when the observed 
degradation signals are complete, fragmented or sparse, and also when the 
scores and errors follow different distributions. Detailed results can be found 
in the supplementary materials [Zhou, Serban and Gebraeel (2010)]. 

7.1.3. Results and analysis of Model 3. We discuss this in the supple- 
mentary materials [Zhou, Serban and Gebraeel (2010)]. 

8. Discussions. In this paper we propose an Empirical Bayesian method 
for predicting the degradation of a partially degraded component or system. 
Specifically, we assume that the degradation process has unknown mean 
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and covariance, which can be estimated through a nonparametric approach 
using a historical database of degradation signals used to estimate the prior 
distribution of the degradation process. These training degradation signals 
may be completely or incompletely observed, that is, may be in the form of 
sparsely observed signals or fragmented signals. 

Our degradation modeling and monitoring approach relies on a series of 
assumptions: 

• The degradation signals follow a Gaussian process. 

• The time points at which the training signals have been observed cover 
the time domain [0,M] cumulatively. 

• The degradation signal of the new component does not cross back the 
threshold D. 

From our simulation results, departures from the Gaussian assumption 
will insignificantly alter the residual life estimates when a large number of 
training signals are observed, as discussed in Section 7. This property is 
inherited from the robustness of the FPCA approach used in estimating the 
empirical prior distribution. 

Under sparse sampling, the selection of the observation times of the train- 
ing degradation signals impacts the accuracy of the degradation prior mod- 
eling. For example, if the degradation signals are uniformly but sparsely 
sampled, the degradation process will not be adequately observed at the 
later extreme time point M, since few components will survive up to this 
time point. Consequently, uniform sampling compromises the accuracy of 
the mean and covariance estimates of the prior degradation process, which, 
in turn, compromises the accuracy of the residual life estimate. In the sim- 
ulation study we show that the accuracy of the residual life estimates is 
low for the traditional uniform sampling as compared to the accuracy of 
the estimates under nonuniform sampling. Thus, the second assumption is 
ensured under nonuniform sampling but not uniform sparse sampling (see 
Section 5). 

The third assumption in our modeling approach relies on that the ex- 
perimenter will shut off or replace the component shortly after it degraded 
beyond the failure threshold D. 

In this paper we have applied the nonparametric approach to crack growth 
data with a wide applicability, for example, in infrastructure (bridges, steel 
structures), maritime (hulls of oil tankers), aeronautical (aircraft fuselage), 
energy (vanes of gas turbines) and others. This case study demonstrates 
the accuracy of the nonparametric approach introduced in this paper as 
compared to random effects parametric models which impose constrains 
on the shape of the trend //(£) and the covariance C(t,t'). Other potential 
applications are relevant to LED data that could be found in Yu and Tseng 
(1998), Liao and Tseng (2006) and Tseng and Peng (2007). 
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SUPPLEMENTARY MATERIAL 

Additional results (DOI: 10.1214/10-AOAS448SUPP; .pdf). In this sup- 
plemental file we provide some additional results of the crack growth data 
study and the simulation study. 
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